#rm(list=ls())
#setwd("/Main Text-Figures/Figure2/")

library(pscl)
library(emIRT)

########################################

alpha <- .5
k <- 0.2
N <- 100
v <- 500
nsim <- 50

full.min <- NULL
full.med <- NULL
full.max <- NULL
sim.info.all <- NULL
rcv.min <- NULL
rcv.med <- NULL
rcv.max <- NULL
full.min.plot <- NULL
full.med.plot <- NULL
full.max.plot <- NULL

rcv.min.plot <- NULL
rcv.med.plot <- NULL
rcv.max.plot <- NULL


for(j in 1:nsim){
  phi <- runif(1,.2,.8)
  
  b1 <- runif(v,0,1)
  b2 <- b1
  sq1 <- runif(v,0,1)
  sq2 <- sq1
  
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  
  d1 <- runif(1,.5,2)
  d2 <- d1
  
########################################
  
  hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
  colnames(hypo.pro.1d) <- c("dim","d","b","sq")
  hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
  colnames(hypo.pro.2d) <- c("dim","d","b","sq")
  hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
  hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
  
########################################
#ideal points
  xL1 <- c(0,runif(npL-1,min=0-d1,max=d1)) #Going to use this that i've put the party leaders in the first and last spots
  xL2 <- xL1
  xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) #Going to use this that i've put the party leaders in the first and last spots
  xR2 <- xR1
  
  left.party <- as.data.frame(cbind(round(xL1,3),round(xL2,3),rep(0,npL)))
  colnames(left.party) <- c("x1","x2","partyid")
  right.party <- as.data.frame(cbind(round(xR1,3),round(xR2,3),rep(1,npR)))
  colnames(right.party) <- c("x1","x2","partyid")
  all.leg <- as.data.frame(rbind(left.party,right.party))
  
########################################
#who votes for it?
  
  vote.func1 <- function(x1,b1,sq1){
    as.numeric(abs(x1-b1)  < abs(x1-sq1) )
  }
  vote.func2 <- function(x2,b2,sq2){
    as.numeric(abs(x2-b2)  < abs(x2-sq2) )
  }
  
  out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
  for(i in 1:N){
    for(j in 1:v){
      out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
    }
  }
  for(i in 1:N){
    for(j in (v+1):(2*v)){
      out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
    }
  }
  
########################################
#does it pass?
  
  passes <- as.numeric(colSums(out) > .5*N)
  hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
  for (i in 1:nrow(hypo.pro.all)){
    hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
  }

########################################
#unity scores  
  left.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
  }
  right.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
  }

  out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] - right.aligned[j]
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]-left.aligned[j]
      } else {
        0
      }
    }
  }

  out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] 
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]
      } else {
        0
      }
    }
  }

########################################
#RCV requested?
  
  out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
  for(i in 1:N){
    for(j in 1:(v*2)){
      out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
    }
  }
  hypo.pro.all$left.cscore <- out.cohesion[1,] 
  hypo.pro.all$right.cscore <- out.cohesion[N,]
  
  hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
  hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
  
  hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
  hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
  
  hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)

########################################
#proposed?
  
  for(i in 1:nrow(hypo.pro.all)){
    hypo.pro.all$proposing.party[i] <- if(hypo.pro.all$b[i] < hypo.pro.all$sq[i]) 0 else 1 
  }
  
  hypo.pro.all$proposed.left <- rep(0,length(nrow(hypo.pro.all)))
  for(i in 1:nrow(hypo.pro.all)){
    if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
    } else {
      0
    }
  }
  hypo.pro.all$proposed.right <- rep(0,length(nrow(hypo.pro.all)))
  for(i in 1:nrow(hypo.pro.all)){
    if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
      hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
    } else {
      0
    }
  }
  hypo.pro.all$proposed <- hypo.pro.all$proposed.left+hypo.pro.all$proposed.right

  all.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1),]
  rcv.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1 & hypo.pro.all$rcv.request==1),]

########################################
  
  vote.matrix.all <- out[,all.proposals$id]
  vote.matrix.rcv <- out[,rcv.proposals$id]

########################################

  y.all <- as.data.frame(vote.matrix.all)
  rc.all <- rollcall(y.all)
  rc.all$votes[rc.all$votes==0]<- -1
  rc.all$codes$nay <- -1
  rcfast.all <- convertRC(rc.all)
  p.all <- makePriors(rcfast.all$n, rcfast.all$m, 1)
  s.all <- getStarts(rcfast.all$n, rcfast.all$m, 1)
  
  fast.out.all <- binIRT(.rc = rcfast.all,
                     .starts = s.all,
                     .priors = p.all,
                     .control = {
                       list(threads = 6,
                            verbose = FALSE,
                            thresh = 1e-6
                       )
                     },
                     .anchor_subject=100
  )
  full.est.all <- round(fast.out.all$means$x,3)
  full.est.all.party <- cbind(seq(1:length(fast.out.all$means$x)),full.est.all,all.leg)
  colnames(full.est.all.party) <- c("id","idp.est.irt","idp1.true","idp2.true","party.id")

  full.min.sim <- length(which(full.est.all.party$idp.est.irt == min(full.est.all.party$idp.est.irt)))
  full.max.sim <- length(which(full.est.all.party$idp.est.irt == max(full.est.all.party$idp.est.irt)))
  full.med.sim <- length(which(full.est.all.party$idp.est.irt == median(full.est.all.party$idp.est.irt))) 
  
  full.min <- c(full.min,as.vector(full.min.sim));full.min
  full.med <- c(full.med,as.vector(full.med.sim),na.rm=TRUE);full.med
  full.max <- c(full.max,as.vector(full.max.sim));full.max

  y.rcv <- as.data.frame(vote.matrix.rcv)
  rc.rcv <- rollcall(y.rcv)
  rc.rcv$votes[rc.rcv$votes==0]<- -1
  rc.rcv$codes$nay <- -1
  rcfast.rcv <- convertRC(rc.rcv)
  p.rcv <- makePriors(rcfast.rcv$n, rcfast.rcv$m, 1)
  s.rcv <- getStarts(rcfast.rcv$n, rcfast.rcv$m, 1)
  
  fast.out.rcv <- binIRT(.rc = rcfast.rcv,
                         .starts = s.rcv,
                         .priors = p.rcv,
                         .control = {
                           list(threads = 6,
                                verbose = FALSE,
                                thresh = 1e-6
                           )
                         },
                         .anchor_subject=100
  )
  full.est.rcv <- round(fast.out.rcv$means$x,3)
  full.est.rcv.party <- cbind(seq(1:length(fast.out.rcv$means$x)),full.est.rcv,all.leg)
  colnames(full.est.rcv.party) <- c("id","idp.est.irt","idp1.true","idp2.true","party.id")

  rcv.min.sim <- length(which(full.est.rcv.party$idp.est.irt == min(full.est.rcv.party$idp.est.irt)))
  rcv.max.sim <- length(which(full.est.rcv.party$idp.est.irt == max(full.est.rcv.party$idp.est.irt)))
  rcv.med.sim <- length(which(full.est.rcv.party$idp.est.irt == median(full.est.rcv.party$idp.est.irt,na.rm=TRUE)))
  
  rcv.min <- c(rcv.min,as.vector(rcv.min.sim));rcv.min
  rcv.med <- c(rcv.med,as.vector(rcv.med.sim));rcv.med
  rcv.max <- c(rcv.max,as.vector(rcv.max.sim));rcv.max
  
  sim.info <- c(nrow(all.proposals),nrow(rcv.proposals),phi,d1)
  sim.info.all <- rbind(sim.info.all,sim.info)
  
}

colnames(sim.info.all) <- c("non.rcv","rcv","phi","d")
sim.info.all <- as.data.frame(sim.info.all)
sim.info.all$percent.rcv <- sim.info.all$rcv / (sim.info.all$rcv+sim.info.all$non.rcv)

pcount <- seq(1,nsim,1)
for(i in 1:nsim){
  fmin <- rep(pcount[i],times=full.min[i])
  fmed <- rep(pcount[i],times=full.med[i])
  fmax <- rep(pcount[i],times=full.max[i])
  rmin <- rep(pcount[i],times=rcv.min[i])
  rmed <- rep(pcount[i],times=rcv.med[i])
  rmax <- rep(pcount[i],times=rcv.max[i])
  full.min.plot <- c(full.min.plot,fmin)
  full.med.plot <- c(full.med.plot,fmed)
  full.max.plot <- c(full.max.plot,fmax)
  rcv.min.plot <- c(rcv.min.plot,rmin)
  rcv.med.plot <- c(rcv.med.plot,rmed)
  rcv.max.plot <- c(rcv.max.plot,rmax)
}


load("idealpoints-fig2-data.RData")

pdf("Fig2_.pdf",width=12,height=8)
par(mfrow=c(2,3),mar=c(3,5,7,1))
hist(full.min.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Minimum Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
mtext(side=2,line=3,"Number of Legislators",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)

hist(full.med.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Median Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)
mtext(side=3,line=4.5,cex=1.75,"Complete Vote Sample")

hist(full.max.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Maximum Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)

par(mar=c(3,5,7,1))
hist(rcv.min.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Minimum Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
mtext(side=2,line=3,"Number of Legislators",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)

hist(rcv.med.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Median Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)
mtext(side=3,line=4.5,cex=1.75,"Roll Call Vote Sample")

hist(rcv.max.plot,breaks=seq(0,nsim,1),col="grey60",border="white",ylim=c(0,100),main="",xlab="",xaxt='n',yaxt='n',ylab="")
mtext(side=3, line=1,"Maximum Ideal Points")
mtext(side=1,line=1,"Simulation #",cex=.75)
mtext(side=2,line=3,"Number of Legislators",cex=.75)
axis(side=1,at=seq(0,50,10),labels=c("1","","","","",50))
axis(side=2,at=c(0,25,50,75,100),labels=c(0,25,50,75,100),las=1)
dev.off()


#median(as.vector(table(rcv.max.plot)))
#median(as.vector(table(full.max.plot)))

#median(as.vector(table(rcv.min.plot)))
#median(as.vector(table(full.min.plot)))
